%% exact sin vs K order consideration
K=100;time=linspace(0,0.5,51);
t1=1:3:199; t2=zeros(size(t1));
%while K<6
exactproj=cell(400,1);
exactproj{1}=@(s) -cos(K*s)/K+1/K;
coeff=zeros(400,1);
for ii=2:400
    a1=K+(ii-1)*pi; a2=(K-(ii-1)*pi);
    exactproj{ii}=@(s) sqrt(2)/2*(1/a1-cos(a1*s)/a1+1/a2-cos(a2*s)/a2);
    coeff(ii)=max(abs(exactproj{ii}(time).^2));
end
hold on
plot((2:400),coeff(2:400));
pause(10);
K=K+3;
%end


%% exact sin vs K order consideration
K=1;time=linspace(0,0.5,51);
t1=1:3:199; t2=zeros(size(t1));
while K<200
exactproj=cell(200,1);
exactproj{1}=@(s) -cos(K*s)/K+1/K;
for ii=2:200
    a1=K+(ii-1)*pi; a2=(K-(ii-1)*pi);
    exactproj{ii}=@(s) sqrt(2)/2*(1/a1-cos(a1*s)/a1+1/a2-cos(a2*s)/a2);
    %exactproj{ii}=@(s) sqrt(2)/(ii-1)/pi*sin((ii-1)*pi*s);
end
termN=1;err=1;
while (err>1/100)
wcevar=zeros(size(time));
for term=1:termN
    wcevar=wcevar+exactproj{term}(time).^2;
end
exactvar=time/2-1/K/4*sin(2*K*time);
err=max(abs(exactvar-wcevar))/max(exactvar);
termN=termN+1;
end
disp(termN-1) 
t2(floor(K/3)+1)=termN-1;
%disp(floor(max(exactvar)*100/pi))
disp(K)
disp(err)
K=K+3;
end
plot(t1,t2)

%% exact sin vs K

K=60;
basemodes=zeros(100,1);
basemodes(2:100)=(1:99)*pi;
amp=sqrt(2);

exactproj=cell(100,1);
exactproj{1}=@(s) -cos(K*s)/K+1/K;
for ii=2:100
    a1=K+(ii-1)*pi; a2=(K-(ii-1)*pi);
    exactproj{ii}=@(s) sqrt(2)/2*(1/a1-cos(a1*s)/a1+1/a2-cos(a2*s)/a2);
    %exactproj{ii}=@(s) sqrt(2)/(ii-1)/pi*sin((ii-1)*pi*s);
end
time=linspace(0,0.5,51); wcevar=zeros(size(time));
for term=1:60
    wcevar=wcevar+exactproj{term}(time).^2;
end
    

hold on
exactvar=time/2-1/K/4*sin(2*K*time);
%exactvar=time;
plot(time,wcevar,'r-');

plot(time,exactvar);





%% dXt=SIGMA(t) W sin(t)
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=2;
err=1; K=5;
while (err>0.01)
% index J for X
alphaK=K; alphas=[zeros(1,alphaK);eye(alphaK)];
fname='alphas'; save(fname,'alphas');
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
% construct orthogoal basis in L2[0,t]
%orthbasis(T); msvalue(num);
% sigma
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=size(alphas,1); X0=zeros(row,1);
% solver of X
disp('solving X')
tic
[time,X] = ode45(@sinvstermrand,tspan,X0,options,[num,alphaK],casenum);
toc

Y=X.*X;
exactvar=time/2-1/4*sin(2*time);
tmp=zeros(size(Y(:,1)));
for ii=1:size(Y,2)
    %tmp=tmp+Y(:,ii).*Y(:,ii);
    tmp=tmp+Y(:,ii);
end

err=max(abs(exactvar-tmp))%/max(exactvar)
K=K+1;
end
K-1

hold on
plot(time,tmp,'r-');
time=linspace(0,0.5,51);
plot(time,exactvar);
title('Variance')
legend('wce var','exact var');

%% dXt=SIGMA(t) W sin(t)
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=3;
err=1; K=10;
while (err>0.01)
% index J for X
alphaK=K; alphas=[zeros(1,alphaK);eye(alphaK)];
fname='alphas'; save(fname,'alphas');
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
% construct orthogoal basis in L2[0,t]
%orthbasis(T); msvalue(num);
% sigma
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=size(alphas,1); X0=zeros(row,1);
% solver of X
disp('solving X')
tic
[time,X] = ode45(@sinvsKrand,tspan,X0,options,[num,alphaK],casenum);
toc

Y=X.*X;
exactvar=time/2-1/8*sin(4*time);
tmp=zeros(size(Y(:,1)));
for ii=1:size(Y,2)
    %tmp=tmp+Y(:,ii).*Y(:,ii);
    tmp=tmp+Y(:,ii);
end

err=max(abs(exactvar-tmp))/max(exactvar)
K=K+1;
end
K-1

hold on
plot(time,tmp,'r-');
time=linspace(0,0.5,51);
plot(time,exactvar);
title('Variance')
legend('wce var','exact var');